A traditional geostatistical workflow to model hydrocarbon reservoirs consists in modeling facies, then populating each facies with petrophysical properties, typically porosity and permeability, using variogram-based algorithms. Because the variogram is a only a two-point measure of spatial variability, variogram-based geostatistics do not allow modeling realistic curvilinear or other geometrically complex facies patterns, such as meandering sand channels, which are critical for connectivity and flow performance assessment.
A more recent modeling approach, referred to as multiple-point statistics simulation, or MPS simulation, has been proposed by Guardiano and Srivastava, Multivariate Geostatistics: Beyond Bivariate Moments: Geostatistics-Troia, in Soares, A., ed., Geostatistics-Troia: Kluwer, Dordrecht, V. 1, p. 133-144, (1993). MPS simulation is a reservoir facies modeling technique that uses conceptual geological models as 3D training images to generate geologically realistic reservoir models. The training images provide a conceptual description of the subsurface geological bodies, based on well log interpretation and general experience in reservoir architecture modeling. MPS simulations extract multiple-point patterns from the training image and anchor the patterns to the well data.
Numerous others publications have been published regarding MPS and its application. Caers, J. and Zhang, T., 2002, Multiple-point Geostatistics: A Quantitative Vehicle for Integrating Geologic Analogs into Multiple Reservoir Models, in Grammer, G. M et al., eds., Integration of Outcrop and Modern Analog Data in Reservoir Models: AAPG Memoir. Strebelle, S., 2000, Sequential Simulation Drawing Structures from Training Images: Doctoral Dissertation, Stanford University. Strebelle, S., 2002, Conditional Simulation of Complex Geological Structures Using Multiple-Point Statistics: Mathematical Geology, V. 34, No. 1. Strebelle, S., Payrazyan, K., and J. Caers, J., 2002, Modeling of a Deepwater Turbidite Reservoir Conditional to Seismic Data Using Multiple-Point Geostatistics, SPE 77425 presented at the 2002 SPE Annual Technical Conference and Exhibition, San Antonio, September 29-October 2. Strebelle, S. and Journel, A, 2001, Reservoir Modeling Using Multiple-Point Statistics: SPE 71324 presented at the 2001 SPE Annual Technical Conference and Exhibition, New Orleans, September 30-October 3.
SNESIM (Single Normal Equation Simulation) is an MPS simulation program which is particularly well known to those skilled in the art of facies and reservoir modeling. In particular, SNESIM simulation is described in detail in Strebelle, S., 2000, Sequential Simulation of Complex Geological Structures Using Multiple-Point Statistics, doctoral thesis, Stanford University and Strebelle, S., 2002, Conditional Simulation of Complex Geological Structures Using Multiple-Point Statistics: Mathematical Geology, V. 34, No. 1. Again, these publications are well-known to facies modelers who employ multiple point statistics in creating facies and reservoir models. These publications of Strebelle are hereby incorporated in their entirety by reference.
Experience has shown that the MPS simulation program SNESIM reproduces training patterns reasonably well. However, SNESIM is significantly more cpu-demanding than a comparable variogram-based simulation program SISIM, also developed at Stanford University. SNESIM requires a very large amount of memory to extract, then store patterns from 3D multimillion-nodes training cubes.
The MPS simulation program SNESIM described in the Strebelle dissertation (2000, pp.40-53) is based on the same sequential simulation paradigm as the traditional indicator variogram-based program SISIM. A condensation of this description SNESIM is contained in Appendix A of this specification. With SISIM the simulation grid nodes are visited one single time along a random path. SISIM is described in Deutsch, C. and Journel, A. (1998) GSLIB: Geostatistical Software Library and User's Guide, second edition, Oxford University Press, New York. Once simulated, a nodal value becomes a hard datum that will condition the simulation of the nodes visited later in the sequence. While in the variogram-based algorithm, kriging is performed at each unsampled node to estimate the local conditional distribution from which a simulated value is drawn, in MPS simulation that local conditional distribution is estimated by scanning a training image with a given data template of conditioning data (Guardiano and Srivastava, 1992).
The main contribution of the Strebelle dissertation (2000) was to decrease the cpu-time of Srivastava's original code by storing ahead of time all required conditional probability distribution functions (cpdf's) in a dynamic data structure called search tree. For convenience, Appendix B of this specification describes how the search tree of Strebelle (2000) is generated. More precisely, denote by W(u) the data search window centered at location u, and τn the data template constituted by the n vectors {hα, α=1 . . . n} defining the n locations u+hα of W(u). Prior to the simulation, the training image is scanned with a data template τn, then the numbers of occurrences of all possible data events associated with data template τn are stored in the search tree. During the MPS simulation, the local cpdf's are retrieved directly from that search tree. Accordingly, the training image need not be scanned anew for each node simulation.
One major limitation of the search tree approach is that the data template τn cannot include too many grid nodes. There are two reasons for such limitation:                1. The amount of memory used to construct the search tree increases exponentially with the size of the data template: for an attribute taking K possible values, e.g. K facies values, the maximum number of possible data events associated with data template τn is Kn. Fortunately that maximum number is rarely reached.        2. The cpu-time needed to retrieve cpdf's from the search tree increases dramatically with a large data template τn. At any unsampled node u, only n′ (≦n) data locations of data template τn are informed with conditioning data (original sample data or previously simulated nodes). Inference of the probability distribution conditional to the data event dn′ constituted by those n′ conditioning data requires calculating the number c(dn′) of occurrences of dn′ in the training image. This number is obtained by summing the numbers of occurrences of all data events dn′ that are associated with data template τn and that include dn′:        
      c    ⁡          (              d        n        ′            )        =            ∑                                                                  d                n                            ⁢                                                          ⁢              associated              ⁢                                                          ⁢              with              ⁢                                                          ⁢                              τ                n                                                                                                        such                ⁢                                                                  ⁢                that                ⁢                                                                  ⁢                                  d                  n                  ′                                            ⋐                              d                n                                                          ⁢          c      ⁡              (                  d          n                )            
The number c(dn) of occurrences of any such data event dn can be read directly from the search tree. The smaller the number n′ of conditioning data, the greater the number of possible data events dn that include dn′, the greater the cpu-time needed to retrieve all the related numbers c(dn′) from the search tree. For an attribute taking K possible values, the number of possible data events dn that include dn′ can be as large as Kn−n′.
The data template cannot include too many grid nodes for memory and cpu-time considerations. Yet, the data search window should be large enough to capture large-scale structures of the training image. One solution to alleviate this conflict is to use a multiple-grid approach, whereby a number G of nested and increasingly finer grids are simulated, see Tran, T., Improving Variogram Reproduction on Dense Simulation Grids, Computers and Geosciences, 20(7):1161-1168 (1994) and Strebelle dissertation (2000, p.46). The g-th (1≦g≦G) grid is constituted by every 2g−1-th node of the final simulation grid. An example of multiple-grid simulation sequence applied to a 2D grid of size 8×8=64 nodes is shown in FIG. 2 with G=3 nested grids.
The data template is re-scaled proportionally to the spacing of the nodes within the grid to be simulated, which allows capturing large-scale training structures with large size templates at the coarse grid levels. The g-th grid is constituted by every 2nd node of the next (g+1)-th grid. In 2D, the number of nodes within the (g+1)-th grid is then 22=4 times the number of nodes within the previously simulated g-th grid. This means that, at the beginning of the simulation of the (g+1)-th grid, ¼ of the nodes within that grid are previously simulated (informed) nodes belonging to the g-th grid. Thus, in order to find at least, say, 20 data to condition the simulation of an unsampled (uninformed) node, a search data template τn constituted by n=4*20=80 nodal locations should be considered. (In 3D, the number of informed nodes in the (g+1)-th grid is 23=8 times the number of uninformed nodes in the g-th grid, and a data template constituted by n=8*20=160 nodes would be needed.) For an attribute taking K=2 possible values, the maximum number of possible data events associated with such template τn is Kn=280=1.3×1024, and the maximum number of possible fully informed events that include a conditioning data event constituted by n′=20 data could be as large as Kn−n′=260=1.2*1018. Fortunately, all these extremely large numbers are capped by the total number of nodes of the training image being scanned.
The present invention addresses previous shortcomings in the computational efficiency of the MPS methodology used in SNESIM.